library(foreign)
library(sandwich)
library(lmtest)
library(ggeffects)

setwd("/Users/joelsievert/Dropbox/Research/Legislative Responsiveness/National Road/data/PRQ Replication")

##########################
##votes to extend road ##
#########################

h21 <- read.csv("h21_extend.csv")
h23 <- read.csv("h23_extend.csv")


#national road
m1 <- glm(rc_natl ~ jackson  + dist_road + dist_natl + slave_p + distrib_cmte, family = binomial, data = subset(h21, slave_state == 1))
o1 <- coefci(m1, vcov = vcovHC(m1, "HC0"), level = 0.9)
o1b <- coefci(m1, vcov = vcovHC(m1, "HC0"), level = 0.95)


#maysville
m2 <- glm(rc_mays ~ jackson + dist_road + dist_mays + slave_p + distrib_cmte, family = binomial, data = subset(h21, slave_state == 1))
o2 <- coefci(m2, vcov = vcovHC(m2, "HC0"), level = 0.9)
o2b <- coefci(m2, vcov = vcovHC(m2, "HC0"), level = 0.95)


##23rd
m3 <- glm(rc ~ jackson + dist_road + dist_natl + slave_p + distrib_cmte , family = binomial, data = subset(h23, slave_state == 1))
o3 <- coefci(m3, vcov = vcovHC(m3, "HC0"), level = 0.9)
o3b <- coefci(m3, vcov = vcovHC(m3, "HC0"), level = 0.95)


####################################
##votes to transfer road to states##
###################################


h17 <- read.csv("h17_transfer.csv")
h20 <- read.csv("h20_transfer.csv")
h23b <- read.csv("h23_transfer.csv")

#17th
m4 <- glm(rc ~  slave_p  + dist + distrib_cmte, family = binomial, data = subset(h17, slave_state == 1))
o4 <- coefci(m4, vcov = vcovHC(m4, "HC0"), level = 0.9)
o4b <- coefci(m4, vcov = vcovHC(m4, "HC0"), level = 0.95)


#20th
m5 <-glm(rc ~ jackson  + dist+ slave_p  + distrib_cmte , family = binomial, data = subset(h20, slave_state == 1))
o5 <- coefci(m5, vcov = vcovHC(m5, "HC0"), level = 0.9)
o5b <- coefci(m5, vcov = vcovHC(m5, "HC0"), level = 0.95)

#23rd

m6 <-glm(rc ~ jackson  + dist+ slave_p  + distrib_cmte , family = binomial, data = subset(h23b, slave_state == 1))
o6 <- coefci(m6, vcov = vcovHC(m6, "HC0"), level = 0.9)
o6b <- coefci(m6, vcov = vcovHC(m6, "HC0"), level = 0.95)


p1 <- ggpredict(m1, terms = "slave_p[0.20, 0.47]", condition = c(jackson = 1, dist_road = 2.1, dist_natl = 4.4, distrib_cmte = 0), ci.lvl = 0.9)
p2 <- ggpredict(m2, terms = "slave_p[0.20, 0.47]", condition = c(jackson = 1, dist_road = 2.1, dist_mays = 5.3, distrib_cmte = 0), ci.lvl = 0.9)
p3 <- ggpredict(m3, terms = "slave_p[0.21, 0.45]", condition = c(jackson = 1, dist_road = 2.7, dist_natl = 6.6, distrib_cmte = 0), ci.lvl = 0.9)
p4 <- ggpredict(m4, terms = "slave_p[0.21, 0.45]", condition = c(jackson = 1, dist = 2.4, distrib_cmte = 0), ci.lvl = 0.9)
p5 <- ggpredict(m5, terms = "slave_p[0.22, 0.48]", condition = c(jackson = 1, dist = 2.2, distrib_cmte = 0), ci.lvl = 0.9)
p6 <- ggpredict(m6, terms = "slave_p[0.2, 0.43]", condition = c(jackson = 1, dist = 2.8, distrib_cmte = 0), ci.lvl = 0.9)




#############
#############
###EXTEND###
############
############

r <- 1:3


par(mfrow = c(1,2))

##extend
par(mar = c(5, 4, 3, 2), bg = "white")
plot(r, c(coef(m1)[5], coef(m2)[5], coef(m3)[5]), ylim = c(-12, 4), xlim = c(0.5, 3.5), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(0.5, 3.5, 0.1), col = 'grey85', lwd = 80)
abline(h = seq(-12, 4, 2), v = seq(0.5, 3.5, 0.5), col = 'white')
abline(h = 0, col = "red")
points(r, c(coef(m1)[5], coef(m2)[5], coef(m3)[5]),  pch = 16)
segments(x1 = r, x0= r, y0 = c(o1[5,1], o2[5,1], o3[5,1]), y1 = c(o1[5,2], o2[5,2], o3[5,2]))

#code puts axis text and an angle
axis(side = 1, at = r, labels = c("Washington", "Maysville", "New Orleans"), lwd = 0 ) 
axis(side = 2, at = seq(-12, 4, 2), las = 2, lwd = 0)

mtext("Estimate", side = 2, line = 2.5)

mtext("Extend Road", side = 3, line = 0.5, cex = 1, font = 2)


#FED CONTROL
par(mar = c(5, 4, 3, 2), bg = "white")
plot(r, c(coef(m4)[2], coef(m5)[4], coef(m6)[4]), ylim = c(-12, 4), xlim = c(0.5, 3.5), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(0.5, 3.5, 0.1), col = 'grey85', lwd = 80)
abline(h = seq(-12, 4, 2), v = seq(0.5, 3.5, 0.5), col = 'white')
abline(h = 0, col = "red")
points(r, c(coef(m4)[2], coef(m5)[4], coef(m6)[4]),  pch = 16)
segments(x1 = r, x0= r, y0 = c(o4[2,1], o5[4,1], o6[4,1]), y1 = c(o4[2,2], o5[4,2], o6[4,2]))

#code puts axis text and an angle
axis(side = 1, at = r, labels = c("17th", "20th", "23rd"), lwd = 0 ) 
axis(side = 2, at = seq(-12, 4, 2), las = 2, lwd = 0)

mtext("Estimate", side = 2, line = 2.5)

mtext("Federal Control", side = 3, line = 0.5, cex = 1, font = 2)

